%%
% 
%  The purpouse of this script is to plot the profiles of the preassure,
%  x velocity and y velocity at 1/8th and 2/8th of the domain after the 
%  expansion.
% 

%%  Load dependencies
addpathsCFD;

%%  Profiles locations

column1 = floor((NIt-NIe)* 1/8);
column2 = floor((NIt-NIe)* 3/8);

%% Plot Structures

[Xp1, Xp2, PT1, PT2, BIP1, BIP2] = PspatialMapPerBlock( h , NIt , NJt , NJe , NIe , PT);
[Xu1, Xu2, UT1, UT2, BIU1, BIU2] = UspatialMapPerBlock( h , NIut , NJut , NJue , NIue , UT);
[Xv1, Xv2, VT1, VT2, BIV1, BIV2] = VspatialMapPerBlock( h , NIvt , NJvt , NJve , NIve , VT); 

%% Plot pressure 

figure(1) 
 
indexes1 = zeros(NJt+2,1);
indexes1(1) = column1 *2-1;
indexes1(end) = column1*2;
for j = 1:NJt 
    indexes1(j+1) = j - 1 + BIP2 + (column1-1) * NJt;
end

indexes2 = zeros(NJt+2,1);
indexes2(1) = column2*2-1;
indexes2(end) = column2*2;
for j = 1:NJt 
    indexes2(j+1) = j - 1 + BIP2 + (column2-1) * NJt;
end

for i=1:floor((T1-T0)/(dt*nSkip))
    title(strcat('pressure at time =',num2str(T0+i*dt*nSkip),'seg.'));
    plot(PT2(indexes1,i),Xp2(indexes1,2),'-b'); hold on;
    plot(PT2(indexes2,i),Xp2(indexes2,2),'-r'); hold off;
    xlabel('Preassure');
    ylabel('Y Axis');
    legend('Profile at 1/8', 'Profile at 3/8');

    maxPT = max([max(max(PT2(indexes1,:))), max(max(PT2(indexes2,:)))]);
    minPT = min([min(min(PT2(indexes1,:))), min(min(PT2(indexes2,:)))]);
    xlim([minPT maxPT]);
    ylim([0 2]);
    pause(0.01);
end

%% Plot velocity x

figure(2) 
 
indexes1 = zeros(NJt+2,1);
indexes1(1) = column1*2-1;
indexes1(end) = column1*2;
for j = 1:NJt 
    indexes1(j+1) = j - 1 + BIU2 + (column1-1) * NJt;
end

indexes2 = zeros(NJt+2,1);
indexes2(1) = column2*2-1;
indexes2(end) = column2*2;
for j = 1:NJt 
    indexes2(j+1) = j - 1 + BIU2 + (column2-1) * NJt;
end

for i=1:floor((T1-T0)/(dt*nSkip))
    title(strcat('Velocity x at time =',num2str(T0+i*dt*nSkip),'seg.'));
    plot(UT2(indexes1,i),Xu2(indexes1,2),'-b'); hold on;
    plot(UT2(indexes2,i),Xu2(indexes2,2),'-r'); hold off;
    xlabel('Velocity x');
    ylabel('Y Axis');
    legend('Profile at 1/8', 'Profile at 3/8');

    maxUT = max([max(max(UT2(indexes1,:))), max(max(UT2(indexes2,:)))]);
    minUT = min([min(min(UT2(indexes1,:))), min(min(UT2(indexes2,:)))]);
    xlim([minUT maxUT]);
    ylim([0 2]);
    pause(0.01);
end

%% Plot velocity y

figure(3) 

indexes1 = zeros(NJt+1,1);
for j = 1:NJt+1 
    indexes1(j) = j-1 + BIV2 + (column1-1) * (NJt+1);
end

indexes2 = zeros(NJt+1,1);
for j = 1:NJt+1 
    indexes2(j) = j-1 + BIV2 + (column2-1) * (NJt+1);
end

for i=1:floor((T1-T0)/(dt*nSkip))
    title(strcat('Velocity x at time =',num2str(T0+i*dt*nSkip),'seg.'));
    plot(VT2(indexes1,i),Xv2(indexes1,2),'-b'); hold on;
    plot(VT2(indexes2,i),Xv2(indexes2,2),'-r'); hold off;
    xlabel('Velocity y');
    ylabel('Y Axis');
    legend('Profile at 1/8', 'Profile at 3/8');

    maxVT = max([max(max(VT2(indexes1,:))), max(max(VT2(indexes2,:)))]);
    minVT = min([min(min(VT2(indexes1,:))), min(min(VT2(indexes2,:)))]);
    xlim([minVT maxVT]);
    ylim([0 2]);
    pause(0.01);
end

%%  Plot all

figure(4) 

indexesP1 = zeros(NJt+2,1);
indexesP1(1) = column1 *2-1;
indexesP1(end) = column1*2;
for j = 1:NJt 
    indexesP1(j+1) = j - 1 + BIP2 + (column1-1) * NJt;
end

indexesP2 = zeros(NJt+2,1);
indexesP2(1) = column2*2-1;
indexesP2(end) = column2*2;
for j = 1:NJt 
    indexesP2(j+1) = j - 1 + BIP2 + (column2-1) * NJt;
end

indexesU1 = zeros(NJt+2,1);
indexesU1(1) = column1*2-1;
indexesU1(end) = column1*2;
for j = 1:NJt 
    indexesU1(j+1) = j - 1 + BIU2 + (column1-1) * NJt;
end

indexesU2 = zeros(NJt+2,1);
indexesU2(1) = column2*2-1;
indexesU2(end) = column2*2;
for j = 1:NJt 
    indexesU2(j+1) = j - 1 + BIU2 + (column2-1) * NJt;
end

indexesV1 = zeros(NJt+1,1);
for j = 1:NJt+1 
    indexesV1(j) = j-1 + BIV2 + (column1-1) * (NJt+1);
end

indexesV2 = zeros(NJt+1,1);
for j = 1:NJt+1 
    indexesV2(j) = j-1 + BIV2 + (column2-1) * (NJt+1);
end

for i=1:floor((T1-T0)/(dt*nSkip))

    subplot(2,2,[2,4]), title(strcat('pressure at time =',num2str(T0+i*dt*nSkip),'seg.'));
    plot(PT2(indexesP1,i),Xp2(indexesP1,2),'-b'); hold on;
    plot(PT2(indexesP2,i),Xp2(indexesP2,2),'-r'); hold off;
    xlabel('Preassure');
    ylabel('Y Axis');
    legend('Profile at 1/8', 'Profile at 3/8');

    maxPT = max([max(max(PT2(indexesP1,:))), max(max(PT2(indexesP2,:)))]);
    minPT = min([min(min(PT2(indexesP1,:))), min(min(PT2(indexesP2,:)))]);
    xlim([minPT maxPT]);
    ylim([0 2]);

    subplot(2,2,1), title(strcat('Velocity x at time =',num2str(T0+i*dt*nSkip),'seg.'));
    plot(UT2(indexesU1,i),Xu2(indexesU1,2),'-b'); hold on;
    plot(UT2(indexesU2,i),Xu2(indexesU2,2),'-r'); hold off;
    xlabel('Velocity x');
    ylabel('Y Axis');
    legend('Profile at 1/8', 'Profile at 3/8');

    maxUT = max([max(max(UT2(indexesU1,:))), max(max(UT2(indexesU2,:)))]);
    minUT = min([min(min(UT2(indexesU1,:))), min(min(UT2(indexesU2,:)))]);
    xlim([minUT maxUT]);
    ylim([0 2]);

    subplot(2,2,3), title(strcat('Velocity x at time =',num2str(T0+i*dt*nSkip),'seg.'));
    plot(VT2(indexesV1,i),Xv2(indexesV1,2),'-b'); hold on;
    plot(VT2(indexesV2,i),Xv2(indexesV2,2),'-r'); hold off;
    xlabel('Velocity y');
    ylabel('Y Axis');
    legend('Profile at 1/8', 'Profile at 3/8');

    maxVT = max([max(max(VT2(indexesV1,:))), max(max(VT2(indexesV2,:)))]);
    minVT = min([min(min(VT2(indexesV1,:))), min(min(VT2(indexesV2,:)))]);
    xlim([minVT maxVT]);
    ylim([0 2]);
    pause(0.01);
end
